% *************************************************************************
% m_randomCBD.m 
% *************************************************************************
clear all; 
close all;
format shortG;
tic;
% *************************************************************************
% ******************** LOAD FUNDAMENTAL PARAMETERS ************************
% *************************************************************************
load parameters.mat;                            
load calib_productivity_amenity.mat;    
load calib_setup.mat;                   
load commuting_cost.mat;  
load calib_auxiliary.mat;

load fundamentals_1950.mat;

load calib_qadjusted.mat; 
load calib_hscale.mat; 

load noagglomeration.mat; 

Rdata=Rdata_all; Ldata=Ldata_all;  
totalpop=sum(Rdata(:,2:end),1);

N=size(ID,1);
NT=size(Rdata,2)-1;           
% *************************************************************************
% ******************** EXOGENOUS FUNDAMENTALS (50-75) *********************
% *************************************************************************
a_mat=[a50, aa];   a_mat=a_mat./repmat(geomean(a_mat),N,1);
b_mat=[b50, bb];   b_mat=b_mat./repmat(geomean(b_mat),N,1);
% *************************************************************************
% ********************** COMMUTING PATTERN 1945 ***************************
% *************************************************************************
Lcom=zeros(N,N,NT+1); 
Lcom(:,:,1)=Lcom45; 
R45=Rdata(:,1); 
L45=Ldata(:,1); 
% *************************************************************************
% ********************** CBD DISTANCE BIN *********************************
% *************************************************************************
cbddist_th=[1; 2; 4; 6]; 
dist45 = zeros(N, length(cbddist_th)+2); 
dist45(:,1) = ID;
dist45(:,2) = cbddist; 
for i=1:N
    if cbddist(i) <= cbddist_th(1);      dist45(i,3) = 1; 
    elseif cbddist(i) > cbddist_th(1) &  cbddist(i) <= cbddist_th(2);  dist45(i,4) = 1;
    elseif cbddist(i) > cbddist_th(2) &  cbddist(i) <= cbddist_th(3); dist45(i,5) = 1;
    elseif cbddist(i) > cbddist_th(3);   dist45(i,6) = 1; 
    end
end
ID_bin1 = dist45(:,1).*dist45(:,3); ID_bin1 = ID_bin1(ID_bin1~=0); n_bin1 = size(ID_bin1,1); 
ID_bin2 = dist45(:,1).*dist45(:,4); ID_bin2 = ID_bin2(ID_bin2~=0); n_bin2 = size(ID_bin2,1);
ID_bin3 = dist45(:,1).*dist45(:,5); ID_bin3 = ID_bin3(ID_bin3~=0); n_bin3 = size(ID_bin3,1);
ID_bin4 = dist45(:,1).*dist45(:,6); ID_bin4 = ID_bin4(ID_bin4~=0); n_bin4 = size(ID_bin4,1);

% *************************************************************************
% *************** MONTECARLO DISTRIBUTION 1945 ****************************
% *************************************************************************
mm_number = 500; 
mm_Rdist = zeros(N,mm_number); 
mm_Ldist = zeros(N,mm_number); 

dist_bin1=zeros(n_bin1,3); dist_bin2=zeros(n_bin2,3); dist_bin3=zeros(n_bin3,3);
for j=1:n_bin1
    dist_bin1(j,:) = [ID_bin1(j), R45(ID_bin1(j)), L45(ID_bin1(j))];
end
for j=1:n_bin2
    dist_bin2(j,:) = [ID_bin2(j), R45(ID_bin2(j)), L45(ID_bin2(j))];
end 
for j=1:n_bin3
    dist_bin3(j,:) = [ID_bin3(j), R45(ID_bin3(j)), L45(ID_bin3(j))];
end 

for m=1:mm_number 

    r_ID_bin4 = ID_bin4(randperm(n_bin4)); dist_bin4=zeros(n_bin4,3);

    for j=1:n_bin4
    dist_bin4(j,:) = [ID_bin4(j), R45(r_ID_bin4(j)), L45(r_ID_bin4(j))];
    end 
    
    dist_bin = [dist_bin1; dist_bin2; dist_bin3; dist_bin4]; 
    dist_bin = sortrows(dist_bin,1);
    mm_Rdist(:,m) = dist_bin(:,2); mm_Ldist(:,m) = dist_bin(:,3); 
end 
%% *************************************************************************
% *************** MONTECARLO DRAW FROM DISTRIBUTION ***********************
% *************** SCALE UP 1950 POPULATION AND EMPLOYMENT *****************
% *************************************************************************
m_R = zeros(N,NT,mm_number); m_L = zeros(N,NT,mm_number); 
for m = 1:mm_number
    m_L(:,:,m) = repmat(mm_Ldist(:,m)./sum(mm_Ldist(:,m)),1,NT).*repmat(totalpop,N,1);
    m_R(:,:,m) = repmat(mm_Rdist(:,m)./sum(mm_Rdist(:,m)),1,NT).*repmat(totalpop,N,1);
end
% ************************************************************************
% ************************** MONTECARLO SET ******************************
% ************************************************************************
mm_success = 0;
maxit=1e+6; maxoutit=1e+4; tol=1e-6; outupdate=0.25; update=0.05; 
mm_Lmat_1950 = []; 
mm_Rmat_1950 = []; 
mm_Rinitial = zeros(N,mm_number); mm_Linitial = zeros(N,mm_number); 

% ************************************************************************
% *************** START MONTECARLO SIMULATION ****************************
% *************************************************************************
for mm=1:mm_number
    
R_i = squeeze(m_R(:,:,mm)); mm_Rinitial(:,mm)=R_i(:,1);  
L_i = squeeze(m_L(:,:,mm)); mm_Linitial(:,mm)=L_i(:,1); 
Q_i = ones(N,NT);

A_i=a_mat.*(L_i./repmat(Area,1,NT)).^alphaest;  
B_i=b_mat.*(R_i./repmat(Area,1,NT)).^betaest;   

Q_n=zeros(N,NT); L_n=zeros(N,NT); R_n=zeros(N,NT); 
H_n=zeros(N,NT+1); H_n(:,1)=hsupply45; H_demand=zeros(N,NT); 
outit=0; 

while outit<maxoutit
    
% continuation values
LQmat=zeros(N,NT); Xmat=zeros(N,NT); Ymat=zeros(N,NT);   
LQmat(:,NT)=Q_i(:,NT); Xmat(:,NT)=B_i(:,NT); Ymat(:,NT)=A_i(:,NT);
    for t=1:NT-1 
    tt=NT-t;
    theta=thetavec(tt); 
    LQmat(:,tt)=Q_i(:,tt).*(LQmat(:,tt+1).^(rho*(1-theta)));
    Xmat(:,tt)=B_i(:,tt).*(Xmat(:,tt+1).^(rho*(1-theta)));
    Ymat(:,tt)=A_i(:,tt).*(Ymat(:,tt+1).^(rho*(1-theta)));
    end 

    for t=1:NT

    if t==1; theta=theta50; Rpre_i=R45; Lpre_i=L45; 
    else     theta=thetavec(t-1); Rpre_i=R_n(:,t-1); Lpre_i=L_n(:,t-1); end 
    
    % conditional probabilities (lambda_L: residential given workplace; lambda_R: workplace given residenace) 
    lambda_L=Kmat(:,:,t).*((repmat(LQmat(:,t)',N,1).^(-mu)).*(repmat(Xmat(:,t)',N,1))).^(rho/sigma); 
    lambda_L=lambda_L./repmat(sum(lambda_L,2),1,N); 
    lambda_R=Kmat(:,:,t).*((repmat(LQmat(:,t),1,N).^(-gammaH/gammaL)).*(repmat(Ymat(:,t),1,N).^(1/gammaL))).^(rho/sigma);
    lambda_R=lambda_R./repmat(sum(lambda_R,1),N,1);
    
    % solve forward equations for population and employment 
    Rpost_i=R_i(:,t); Lpost_i=L_i(:,t);
    it=0;
    while it<maxit 
    Rpost_n=(1-theta).*Rpre_i+lambda_L'*(Lpost_i-(1-theta).*Lpre_i);
    Lpost_n=(1-theta).*Lpre_i+lambda_R*(Rpost_i-(1-theta).*Rpre_i);  
    if max(abs(Rpost_n-Rpost_i)) < tol && max(abs(Lpost_n-Lpost_i)) < tol; break;
    else
       Rpost_i=update.*Rpost_n+(1-update).*Rpost_i; 
       Lpost_i=update.*Lpost_n+(1-update).*Lpost_i;
       it=it+1;  
    end
    end
    R_n(:,t)=Rpost_n; L_n(:,t)=Lpost_n;

    % commuting patterns (column:=workplace; row:=residential place)
    Lcompost=(1-theta).*Lcom(:,:,t)+lambda_R.*repmat((Rpost_n-(1-theta).*Rpre_i)',N,1); 
    Lcom(:,:,t+1)=Lcompost;
    
    % wage (workplace) and average income (residential place) 
    w_n=(Q_i(:,t).^(-gammaH/gammaL)).*(A_i(:,t).^(1/gammaL));   
    v_n=(Lcompost./repmat(Rpost_n',N,1))'*w_n; 
    
    % aggregate demand and supply (stock) of floor space
    floor_demand=mu.*v_n.*Rpost_n + (gammaH/gammaL).*w_n.*Lpost_n;
    h_n=(1-psi).*H_n(:,t)+hbar(t)*(Q_i(:,t).^eta).*Area;
    H_n(:,t+1)=h_n; 
    H_demand(:,t)=floor_demand; 

    % floor space market clearing condition
    Q_n(:,t)=H_demand(:,t)./H_n(:,t+1); 
    end

A_n=a_mat.*(L_n./repmat(Area,1,NT)).^alphaest;  
B_n=b_mat.*(R_n./repmat(Area,1,NT)).^betaest;   
dd_Q=Q_n-Q_i; dd_A=A_n-A_i; dd_B=B_n-B_i; 


if mod(outit,20)==1
   disp(['Montecarlo is:', ' ', num2str(mm), ' ', 'and Largest Error is:', ' ',num2str(outit), ' ',num2str(max(abs(dd_Q)))]);
end 

if max(max(abs(dd_A)))<tol && max(max(abs(dd_B)))<tol && max(max(abs(dd_Q)))<tol 
   mm_Lmat_1950 = [mm_Lmat_1950, L_n(:,1)]; 
   mm_Rmat_1950 = [mm_Rmat_1950, R_n(:,1)];
   mm_success = mm_success+1; 
   break
elseif outit == maxoutit
   break
else   
   outit = outit+1;
   Q_i=update.*Q_n+(1-update).*Q_i;
   B_i=update.*B_n+(1-update).*B_i; 
   A_i=update.*A_n+(1-update).*A_i; 
end

end
end

%% *************************************************************************
% ******************* MONTECARLO RESULTS **********************************
% *************************************************************************
ldens_R_mm = log(mm_Rmat_1950./repmat(Area,1,mm_success)); 
ldens_L_mm = log(mm_Lmat_1950./repmat(Area,1,mm_success)); 
[Max_Rdens, Max_Rdens_ID] = max(ldens_R_mm);
[Max_Ldens, Max_Ldens_ID] = max(ldens_L_mm);
MC_result_R=zeros(N,1); 
MC_result_L=zeros(N,1); 
for i=1:mm_success
    MC_result_R = MC_result_R + (Max_Rdens_ID(i)==ID); 
    MC_result_L = MC_result_L + (Max_Ldens_ID(i)==ID); 
end
MCprob_result_R = MC_result_R./mm_success; 
MCprob_result_L = MC_result_L./mm_success; 


save('rationalpath_randomCBD.mat', 'mm_Rmat_1950', 'mm_Lmat_1950');

result=array2table([ID, cbddist, Area, R45,L45, Rdata(:,2), Ldata(:,2), MC_result_R, MC_result_L,MCprob_result_R,MCprob_result_L],'VariableNames',{'geocode','cbddist', 'area',...
     'pop45','emp45','pop50','emp50',...
      'MC_R','MC_L','MCprob_R','MCprob_L'}); 
writetable(result,'../../output/quant/output_randomcbd.csv');

disp(['Elapsed time is: ', num2str(toc / 60), ' minutes']);